# Data from USGS National Water Information System (NWIS)

# Arizona state shape (CRS = 5070)
az = us_states() %>% 
  filter(name == 'Arizona') %>% 
  st_transform(5070)

# add 'well #' for each unique well
usgs = az_nwis_unique_sites
usgs$wellid = paste("well", 1:nrow(usgs))

# remove wells located at same lat/long but have different well ID
usgs_spatial = usgs %>% 
  st_as_sf(coords = c('long_nad83', 'lat_nad83'), crs = 4269) %>% 
  st_transform(5070) %>%
  as_Spatial() %>% 
  remove.duplicates() %>% 
  st_as_sf() %>% 
  group_by(wellid)

# match 'well #' from spatial data frame to corresponding well in time series data frame
usgs_time = az_nwis_all %>% 
  group_by(site_id) 
usgs_time = left_join(usgs_time, select(usgs, site_id, wellid), by = "site_id") 
usgs_time = usgs_time %>% filter(wellid %in% usgs_spatial$wellid)

# FILTER BY THRESHOLD - at least 10 measurements and at least 1 measurement every 5 years
thresh = usgs_time %>% 
  group_by(wellid) %>% 
  filter(measurement_dist >= 10) %>% 
  mutate(measure_period = year - lag(year)) %>% 
  filter(any(measure_period > 5))
usgs_time = usgs_time %>% filter(!wellid %in% thresh$wellid, measurement_dist >= 10) %>% 
  mutate(measure_period = year - lag(year))
usgs_spatial = usgs_spatial %>% 
  filter(!wellid %in% thresh$wellid) %>% 
  filter(measurement_dist >= 10)


# Data from Arizona Department of Water Resources (ADWR)

# add 'well #' for each unique well
state = adwr_unique_sites %>%
  filter(date > 1960-01-01, date < as.Date.character('2000-01-01'))
state$wellid = paste("well", 1:nrow(state))

# remove wells located at same lat/long but have different well ID
state_spatial = state %>%
  st_as_sf(coords = c('long_nad83', 'lat_nad83'), crs = 4269) %>%
  st_transform(5070) %>%
  as_Spatial() %>%
  remove.duplicates() %>%
  st_as_sf() %>%
  group_by(site_id)

# match 'well #' from spatial data frame to corresponding well in time series data frame
state_time = adwr_all %>%
  group_by(site_id)
state_time = left_join(state_time, select(state, site_id, wellid), by = "site_id")
state_time = state_time %>% filter(wellid %in% state_spatial$wellid)

# FILTER BY THRESHOLD - at least 10 measurements and at least 1 measurement every 5 years
thresh_state = state_time %>%
  group_by(wellid) %>%
  filter(measurement_dist >= 10) %>%
  mutate(measure_period = year - lag(year)) %>%
  filter(any(measure_period > 5))
state_time = state_time %>% filter(!wellid %in% thresh_state$wellid, measurement_dist >= 10) %>%
  mutate(measure_period = year - lag(year))
state_spatial = state_spatial %>%
  filter(!wellid %in% thresh_state$wellid) %>%
  filter(measurement_dist >= 10)

# PLOTS A BUFFER AROUND SPECIFIED WELL AND LOCATES OTHER WELLS INSIDE THE BUFFER AREA
buffer_fun = function(df, well, buff, state) {
  buffer = st_buffer(df[well,], buff)
  near1 = st_intersection(df[,], buffer) %>% filter(measurement_dist >= 10)
  
  plot = ggplot() + 
    geom_sf(data = state) +
    geom_sf(data = buffer, fill = NA) + 
    geom_sf(data = near1, col = "red", size = .5) + 
    labs(caption = paste(nrow(near1), 'wells')) +
    theme_void() +
    theme(plot.caption = element_text(size = 22, face = "bold", hjust = 0.5))
  print(plot)
  return(near1)
}

# RETURNS DATAFRAME OF WELLS IN DTW RANGE (MIN - MAX)
dtw_range = function(df, min, max) {
  df = df %>% filter(dtw <= max, dtw >= min) %>% 
    mutate(sd = sd(dtw)) %>% 
    arrange(desc(date))
  return(df)
}


# PLOTS AN INDIVIDUAL WELL BY WELL ID (DTW TIME SERIES)
plot_well = function(df_time, id) {
  wells = df_time %>% filter(wellid == paste('well', id))
  plot = ggplot(data = wells, aes(x = date, y = dtw_ft)) +
    geom_line(aes(y = dtw_ft, col = wellid), size = 2) +
    scale_y_reverse() +
    labs(title = paste('Well', id),
         caption = paste('Min depth:', min(wells$dtw_ft), '\nMax depth:',
                          max(wells$dtw_ft), '\nNumber of measurements:',
                          wells$measurement_dist),
         x = 'Year',
         y = 'DTW (ft)') + 
    theme_bw() +
    theme(plot.title = element_text(face = 'bold',color = 'black', size = 18, hjust = 0.5),
          axis.text.x = element_text(color="black", size=14), 
          axis.text.y = element_text(color="black", size=14), 
          axis.title.x = element_text(face = 'bold', color="black", size=16), 
          axis.title.y = element_text(face = 'bold', color="black", size=16),
          plot.caption = element_text(face = 'bold', color = 'black', size = 14),
          panel.grid.major = element_line(colour = "#808080"),
          panel.grid.minor = element_line(colour = "#808080", size = 1),
          legend.position = 'none') 
  print(plot)
}

# PLOTS MULTIPLE WELLS IN DTW RANGE
multi_well_plot = function(df, df_time, min, max) {
  df = df %>% filter(dtw <= max, dtw >= min) %>% 
    arrange(desc(date))
  df_time = df_time %>% filter(wellid %in% df$wellid)
  plot1 = ggplot(data = df_time, aes(x = date, y = dtw_ft)) +
    geom_line(aes(y = dtw_ft, col = wellid), size = 1) +
    scale_y_reverse() +
    labs(x = 'Year',
         y = 'DTW (ft)') + 
    theme_bw() +
    theme(plot.title = element_text(face = 'bold',color = 'black', size = 18, hjust = 0.5), 
          axis.text.x = element_text(color="black", size=14), 
          axis.text.y = element_text(color="black", size=14), 
          axis.title.x = element_text(face="bold", color="black", size=16), 
          axis.title.y = element_text(face="bold", color="black", size=16), 
          panel.grid.major = element_line(colour = "#808080"),
          panel.grid.minor = element_line(colour = "#808080", size = 1),
          legend.position = 'none')
  print(plot1)
  return(df)
}


# PLOTS WELLS DTW WITHIN A SPECIFIED RANGE, COLOR GRADIENT BY DEPTH, RETURNS CORRESPONDING DATA FRAME
map_dtw = function(state, df, min, max) {
  df = df %>% filter(dtw <= max, dtw >= min) %>% 
    mutate(sd = sd(dtw)) %>% 
    arrange(dtw)
  plot1 = ggplot() + 
    geom_sf(data = state, size = 1, col = 'black') +
    geom_sf(data = df, aes(col = dtw), size = 1.5) + 
    theme_void() +
    theme(plot.title = element_text(face = "bold", color = "black", size = 18),
          plot.subtitle = element_text(size = 12),
          plot.caption = element_text(face = 'bold', size = 12), 
          legend.title = element_text(color="black", size=14), 
          legend.text = element_text(color="black", size=14)) +
    scale_colour_gradient() +
    labs(title = paste('DTW', min, '-', max, 'ft'),
         fill = 'Number of wells', 
         col = 'Depth to Water (ft)')
  print(plot1)
  return(df)
}